1. Field of the Invention
The present invention relates generally to seismic signal processing and, more particularly, to a method for efficient inversion of near singular geophysical signals.
2. Description of the Related Art
A common data processing activity in exploration geophysics is inverting geophysical signals. This can arise in Wiener filtering in its various forms, for example, deconvolution, cross equalization, adaptive filtering and echo cancellation. Such signals are characterized by a convolutional equation system and inverting the signals are all variations on the convolutional model least squares problem. These problems are discussed in Treitel, S. and Lines, L. R., xe2x80x9cLinear Inverse Theory and Deconvolutionxe2x80x9d Geophysics, vol. 47, p p. 1153-1159, 1982 and in Lines, L. R. and Treitel, S., xe2x80x9cTutorial: A Review of Least-Squares Inversion and its Applications to Geophysical Problemsxe2x80x9d, Geophysical Prospecting, Vol. 32, pp.159-186, 1984. For geophysical signals, the corresponding systems are often near singular (ill-conditioned) or rank-deficient. Classic methods for inverting such systems require stabilization methods such as adding white random noise, using a xcex3 damping factor, spectral shaping or constraint, or other so-called regularization methods. For the convolutional model, this is usually done in the normal equation form using Levinson recursion on the Toeplitz autocorrelation matrix, as will be explained below.
In this class of problems, one is given an input signal i(t), an output signal o(t) and a finite impulse response filter u(t) having the relationship:
i(t)=i {circle around (x)}u(t)=o(t)xe2x80x83xe2x80x83(1)
Here the symbol {circle around (x)} indicates convolution and the vectors i, u and o are each functions of time t. It is typically desired to find the particular solution filter or wavelet û(t) which which gives the minimum error e(t) in the least squares sense. The solution to the least squares minimization problem is found by first defining the error vector by
e(t)=i(t){circle around (x)}ûl (t)xe2x88x92o(t)xe2x80x83xe2x80x83(2)
and then minimizing the mean squared scalar error xcex5 defined by the inner product of the error vector e from Eq. (2) with itself
xcex5=eHxc2x7e.xe2x80x83xe2x80x83(3)
Here the superscript H indicates the Hermetian transpose operator or complex conjugate, and the dot xc2x7 indicates the vector inner product.
For a filter vector u(t) of length n and an output vector o(t) of length m in Eq. (1), the input vector i(t) will be of length m+nxe2x88x921. Convolution of the input vector i with filter vector u may be represented by multiplication of u by an mxc3x97n matrix T, where T is the Toeplitz matrix defined by
T(t,xcfx84)=i(nxe2x88x92xcfx84+txe2x88x921)xe2x80x83xe2x80x83(4)
for t=0, . . . m+nxe2x88x922 and xcfx84=0, . . . nxe2x88x921. FIG. 1 shows the matrix T. Eq. (4) allows Eq. (1) to be rewritten in matrix form as
Txc2x7u=oxe2x80x83xe2x80x83(5)
Eq. (5) is known as the direct form equations of the least squares problem. In geophysical signal processing, the number of data observations, represented by rows m, is usually larger than the number of physical parameters, represented by columns n. Thus the direct form equations, Eq. (5), are typically a system of m overdetermined equations.
The least squares solution û may be found from the normal form equations, an exactly determined Toeplitz nxc3x97n system
(THxc2x7T)xc2x7û=THxc2x7oxe2x80x83xe2x80x83(6)
Eq.(6) is Hermetian Toeplitz if i(j)=0 for j=0 . . . , nxe2x88x922 and j=m, . . . , m+nxe2x88x922. This is shown in FIG. 2, which indicates a column circular matrix. FIG. 1 then represents the transient-free case, while FIG. 2 represents the more typical fully-windowed case.
Eq. (6) is typically simply written as
Axc2x7û=cxe2x80x83xe2x80x83(7)
where A=THxc2x7T is the nxc3x97n autocorrelation or Gram matrix and c=THxc2x7o is the nxc3x971 cross-correlation vector. The typical method for solving Eq. (7) has been Levinson recursion, which will solve an nxc3x97n Toeplitz system proportional to n2 multiply-add steps, that is, has order O(n2) computational complexity. The computations of matrix A and vector c will then add another (m+1)xe2x88x92n steps to these calculations.
The convolutional model least squares problems encountered in geophysical signal processing are often ill-conditioned. A measure of the condition of a matrix T is the condition number cond(T). The condition number of a matrix T is defined as the ratio of the largest to the smallest singular values "sgr"(i) in the singular value decomposition of T or the square root of the ratio of the largest to the smallest eigenvalues xcex(i) in the eigenvalue decomposition of the corresponding Gram matrix A. Thus
cond(THxc2x7T)=cond(T)2=cond(A)=("sgr"max/"sgr"min)2=xcexmasx/xcexminxe2x80x83xe2x80x83(8)
where
"sgr"max=maximum singular value of T,
"sgr"min=minimum singular value of T,
xcexmax=maximum eigenvalue of A, and
xcexmin=minimum eigenvalue of A.
Thus the normal form equations given in Eqs. (6) and (7) have a condition number which is the square of the condition number of the direct form equations of Eq. (5).
Matrices with larger condition numbers are more difficult to accurately invert. This suggests that it would be advantageous to solve the direct form of the least squares problem in Eq. (5) without resorting to conversion to the less well-conditioned normal form of Eqs. (6) and (7). This is usually done by using any number of matrix decomposition techniques, all of which are O(mxc2x7n2) in computational complexity. The advantage to the normal form equation approach for Toeplitz systems is reduced computations at the expense of numerical stability.
The matrix A is defined as singular if xcexmin=0, and as nearly singular or ill-conditioned if |xcexmin| less than xcex4, for an appropriately selected small xcex4. The eigenvalues of Gram matrices are nonnegative. For a practical lower limit, one may choose xcex4=nxc2x7xcexmaxxc2x7eps where eps is the limit of relative floating point accuracy on the computer used. For IEEE standard 754 double precision floating point arithmetic, eps=2xe2x88x9252≈10xe2x88x9216 The maximum conditioning number defining ill-conditioning, cmax, may be estimated by substituting |xcexmin| less than xcex4 into Eq. (8), yielding cmax=1/(nxc2x7eps). Then a nearly singular or ill-conditioned system is defined by cond(A) greater than cmax.
For a well-conditioned matrix A, a least squares solution û to Eq. (7) is typically found with a matrix decomposition technique such as QR or LU matrix decomposition. The ill-conditioned matrix case, where cond(A) greater than cmax, may be solved using singular value decomposition. The rank-deficient case may also be solved by the singular value decomposition approach. Another approach for an ill-conditioned matrix A is to use a small damping factor satisfying xcex3 less than  less than |xcexmax| and xcex4 less than xcex3 less than |xcexmin|. This leads to solving the new system
(A+xcex3xc2x7I)xc2x7û=cxe2x80x83xe2x80x83(9)
where I is the identity matrix. The singular value decomposition is the stablest approach, but also the slowest approach.
In geophysical applications, the spectral content of the input signal i(t) is often characterized by regions of low amplitude caused by natural shaping, filtering, and noise or reverberation. Such spectral gaps or notches make spectral inversion very difficult without resorting to stabilization methods such as adding white random noise, using a xcex3 damping factor as in Eq.(9), spectral shaping or constraint, or other so-called regularization methods. The solutions, wavelets, obtained by traditional methods are often ill-behaved and as well may not give useful results. Whereas using methods such as singular value decomposition lead to excellent results, they are generally too computationally expensive to be used in most situations.
U.S. Pat. No. 5,864,786 to Jericevic and having the same assignee as the present invention, describes a method of spectral pruning to solve linear systems representing physical systems. However, that application does not take advantage of using the corresponding Toeplitz matrix or solving the corresponding normal form equations to deal with ill-conditioned systems.
The present invention is a method for inverting a near singular geophysical signal. First, an input seismic signal is represented as an input vector i(t) of length m+nxe2x88x921, an output seismic signal is represented as an output vector o(t) of length m, and a finite impulse response filter is represented as a solution vector u(t) of length n, satisfying a convolutional equation i(t){circle around (x)}u(t)=o(t). Next an mxc3x97n Toeplitz matrix T (t,xcfx84) corresponding to i(t) is calculated, satisfying a matrix equation T (t,xcfx84)xc2x7u(t)=o(t). Both sides of the matrix equation are transformed to the frequency domain, to generate a transformed equation. Then, both sides of the transformed equation are spectrally pruned, to generate a pruned equation. Finally the pruned equation is solved.